function Fig4_a_b_c_d()
lb = 25.66 / 6^0.5;
EC = 56.1 * 6^0.5;

% These files are created by process_figure4_data.ipynb
% process_figure4_data.ipynb is located in .../processed data/
data_A_r = readmatrix('processed data/multi_filling_lineC_r.csv');
r_A = data_A_r(1,4:end)*lb;
nu_A = data_A_r(2:end,3);
phi_r_A = data_A_r(2:end,4:end)'*EC;

% These files are created by process_figure4_data.ipynb
% process_figure4_data.ipynb is located in .../processed data/
data_A_q = readmatrix('processed data/multi_filling_lineC_q.csv');
q_A = data_A_q(1,4:end);
phi_q_A = data_A_q(2:end,4:end)';

% These files are created by process_figure4_data.ipynb
% process_figure4_data.ipynb is located in .../processed data/
data_B_r = readmatrix('processed data/multi_filling_lineB_r.csv');
r_B = data_B_r(1,4:end)*lb;
nu_B = data_B_r(2:end,3);
phi_r_B = data_B_r(2:end,4:end)'*EC;

% These files are created by process_figure4_data.ipynb
% process_figure4_data.ipynb is located in .../processed data/
data_B_q = readmatrix('processed data/multi_filling_lineB_q.csv');
q_B = data_B_q(1,4:end);
phi_q_B = data_B_q(2:end,4:end)';


Fcolor='k';
Bcolor='w';
Fsize=18;  %font size

f = figure();
f.Color = Bcolor;
%% plot
% f = figure()
% f.Color = Bcolor;

subplot(2,2,1)
imagesc(nu_A,r_A,phi_r_A);

ax=gca;
ax.Color=Bcolor;
ax.XColor=Fcolor;ax.YColor=Fcolor;
xl=xlabel('filling factor');
yl=ylabel('r (nm)');
ax.FontSize=Fsize;
yl.FontSize=Fsize;xl.FontSize=Fsize;
ax.YDir='normal';ax.XDir='normal';

cMap=getColor('Magma');colormap(ax,cMap);
ax.CLim=[-0.8,0.4];
axis([-0.8,-0.2,0,max(r_A)])
%%
% f = figure()
% f.Color = Bcolor;
subplot(2,2,2)
imagesc(nu_A,q_A,phi_q_A);

ax=gca;
ax.Color=Bcolor;
ax.XColor=Fcolor;ax.YColor=Fcolor;
xl=xlabel('filling factor');
yl=ylabel('k lB');
ax.FontSize=Fsize;
yl.FontSize=Fsize;xl.FontSize=Fsize;
ax.YDir='normal';ax.XDir='normal';

cMap=getColor('viridis');colormap(ax,cMap);
ax.CLim=[-0.12,0.02];
axis([-0.8,-0.2,0.5,3])
%%
% f = figure()
% f.Color = Bcolor;
subplot(2,2,3)
imagesc(nu_B,r_B,phi_r_B);

ax=gca;
ax.Color=Bcolor;
ax.XColor=Fcolor;ax.YColor=Fcolor;
xl=xlabel('filling factor');
yl=ylabel('r (nm)');
ax.FontSize=Fsize;
yl.FontSize=Fsize;xl.FontSize=Fsize;
ax.YDir='normal';ax.XDir='normal';

cMap=getColor('Magma');colormap(ax,cMap);
ax.CLim=[-1,0.5];
axis([-0.8,-0.2,0,max(r_A)])
%%
% f = figure()
% f.Color = Bcolor;
subplot(2,2,4)
imagesc(nu_B,q_B,phi_q_B);

ax=gca;
ax.Color=Bcolor;
ax.XColor=Fcolor;ax.YColor=Fcolor;
xl=xlabel('filling factor');
yl=ylabel('k lB');
ax.FontSize=Fsize;
yl.FontSize=Fsize;xl.FontSize=Fsize;
ax.YDir='normal';ax.XDir='normal';

cMap=getColor('viridis');colormap(ax,cMap);
ax.CLim=[-0.25,0.04];ax.CLim=[-0.12,0.02];
axis([-0.8,-0.2,0.5,3])
end
